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Results presented at the Fifth Drag Prediction Workshop using CFL3D, FUN3D, and 
NSU3D are described. These are calculations on the workshop provided grids and drag 
adapted grids. The NSU3D results have been updated to reflect an improvement to skin 
friction calculation on skewed grids. FUN3D results generated after the workshop are 
included for custom participant generated grids and a grid from a previous workshop. Uni- 
form grid refinement at the design condition shows a tight grouping in calculated drag, 
where the variation in the pressure component of drag is larger than the skin friction com- 
ponent. At this design condition, A fine-grid drag value was predicted with a smaller drag 
adjoint adapted grid via tetrahedral adaption to a metric and mixed-element subdivision. 

The buffet study produced larger variation than the design case, which is attributed to 
large differences in the predicted side-of-body separation extent. Various modeling and 
discretization approaches had a strong impact on predicted side-of-body separation. This 
large wing root separation bubble was not observed in wind tunnel tests indicating that 
more work is necessary in modeling wing root juncture flows to predict experiments. 

I. Introduction 

The AIAA Drag Prediction Workshop (DPW) series has been an invaluable forum to evaluate Computa- 
tional Fluid Dynamics (CFD) tools. A detailed history of the DPW series and its objective are available from 
Vassberg et al. 1 and Morrison. 2 An abbreviated history and summary of key results are provided here to 
provide context for this DPW-V contribution. The objectives for these workshops include the assessment of 
the the state-of-the-art computational methods as practical aerodynamic tools for aircraft force and moment 
prediction of industry relevant geometries, the provision of an impartial forum for evaluating the effectiveness 
of existing computer codes and modeling techniques using Reynolds averaged Navier-Stokes (RANS) solvers, 
and identification of areas needing additional research and development. 

DPW-I 3 solicited CFD predictions of the lift, drag, and pitching moment for a subsonic transport wing 
body, the DLR.-F4. The conclusions of Levy et al. 3 set the tone for the subsequent workshops, “The overall 
level of scatter is too high and needs to be reduced to determine overall accuracy and trends with grid type, 
turbulence model, etc. Future work should try to identify sources of the scatter, for example, grid quality.” 
DPW-II 4 focused on grid refinement studies and the prediction of installed pylon-nacelle drag increments for 
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the DLR-F6 wing-body and wing-body-nacelle-pylon configurations. The DLR-F6 has substantial areas of 
separation at the wing-body juncture and at the wing-pylon juncture. These regions of flow separation were 
considered a likely culprit for the observed lack of grid convergence. 4,5 DPW-III 6 introduced two wing only 
configurations and retained the DLR-F6 wing-body with the addition of a new wing-body fairing designed 
to eliminate the side-of-body separation. Details of the FX2B fairing design are documented by Vassberg, 
Sclafani, and DeHann' and verified in a wind tunnel test. 8 A higher Reynolds number was specified in 
DPW-III to reduce the wing upper-surface trailing-edge flow separation. DPW-III was held between the 
design and test of the fairing, so this workshop was truly a prediction. After the workshop was complete, 
the test data was compared 9 to the results submitted to DPW-III. 

The DPW-III cases showed a reduction in variation with grid refinement. 10 However, the configuration 
including the FX2B fairing without wing-root separation had a similar level of variation as the the DLR- 
F6 configuration with wing-root separation. Eliminating the wing root separation region did not reduce 
variation. Morrison and Hemsch 10 conclude, “After three workshops, it is still clear that grids remain a 
first order effect and obtaining high quality grids is the first step to obtaining a high quality solution. 
Furthermore, obtaining a high quality family of grids for relatively simple wing-body configurations suitable 
for grid convergence studies remains a formidable challenge.” Grids were considered a first order effect even 
with the focused attempt to make a series of grids with uniform refinement to grids that were much larger 
than those used in previous workshops. A reoccurring theme of grid quality and resolution issues during the 
first three workshops was identified and documented by Mavriplis et al. 11 

II. Common Research Model 

For DPW-IV, a new configuration, the Common Research Model (CRM) was created. This wing body 
configuration with optional horizontal tail and nacelle with pylons is designed to be representative of a 
contemporary high performance transonic transport. 12 The CRM without nacelles was analyzed with and 
without a horizontal tail in DPW-IV. 1 The CRM without nacelles or horizontal tail is the focus of DPW- 
V. The derived reference quantities of the full-scale vehicle are summarized in Table 1, which correspond 
to the geometry and grids provided by the DPW committee. A wind tunnel test was performed on this 
configuration at the NASA National Transonic Facility (NTF). 13 Additional tests were performed at the 
NASA Ames 11-ft wind tunnel 14 to quantify the repeatability of the measurements. The NTF data is 
denoted as NTF-R.UN-44 and the Ames 11-ft data is denoted 11FT-R.UN-126 in Fig. 1 and Fig. 2. The 
NASA Ames 11-ft wind tunnel tests included pressure sensitive paint 15 and skin friction measurements. 16 
The model support system is not included in the DPW-IV or DPW-V computational model. Including the 
support system causes a shift in pitching moment as noted in a computational investigation. 1 ' This was 
verified in an additional computational investigation, 18 where adjusting the wing twist caused an additional 
shift in predicted pitching moment toward the wind tunnel measurements. 


Table 1. Reference geometry for the CRM. 


Mean Aerodynamic Chord ( c re f ) 

275.80 inch 

Wing Reference Area/2 

297, 360 inch 2 

Wing Span/ 2 

1,159.75 inch 

X Moment Center 

1,325.9 inch 

Z Moment Center 

177.95 inch 

Aspect Ratio(AR) 

9.0 


III. DPW-IV 

The DPW-IV summary by Vassberg et al. 1 shows a large increase in code-to-code variation in lift and 
pitching moment between 3 and 4 degrees angle of attack. This large increase in variation was attributed 
to a large range of predicted side-of-body separation size. This DPW-V study also indicates that differences 
in predicted side-of-body separation extent at angles of attack above the design condition dramatically 
increased the scatter of results. At DPW-IV, the participants were unsure of the validity of this side-of-body 
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(a) Coefficient of lift as a function of angle of attack. (b) Coefficient of drag as a function of angle of attack. 

Figure 1. Measured CRM wing-body forces without nacelles or horizontal tail. 



Figure 2. Measured CRM wing-body lift Coefficient as a function of pitching moment coefficient without nacelles or 
horizontal tail. 
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separation because the configuration was not wind tunnel tested before DPW-IV. The wind tunnel tests 13, 14 
after DPW-IV did not indicate the abrupt drop in lift associated with the side-of-body separation predicted 
by some of the DPW-IV participants. Neither measured skin friction or oil flow visualizations 16 indicated 
a large wing root separation bubble at 4 degrees angle of attack. Pressure-sensitive paint measurements 15 
indicated the pressure distributions predicted by the computational methods without massive a separation 
bubble. DPW-IV included a horizontal tail for most cases, except one tail-off case for the downwash study. 
Sclafani et al. 19 showed that the presence of the tail or varying the tail incidence has a smaller influence on 
side-of-body separation than other factors. 

The DPW-IV participants and others published evidence of a large range of predicted side-of-body 
separation sizes. Sclafani et al. 19 exercised a number of OVERFLOW discretization options on different grid 
topologies. A bubble was shown at 4 degrees angle of attack for an upwind scheme but not a central scheme. 
Modeling the Full Navier-Stokes (FNS) viscous terms had a larger bubble than a Thin-Layer Navier-Stokes 
(TLNS) approximation for a number of different grids at 4 degrees angle of attack. Increasing the wing root 
juncture grid resolution with a finer orthogonal grid increased the extent of separation. Mani et al. 20 used 
OVERFLOW to show a small bubble at the CRM design condition that shrunk with grid refinement. The 
bubble was larger with FNS than the TLNS approximation. BCFD simulations on a grid with more wing 
root juncture refinement showed a larger bubble with TLNS that also shrank with grid refinement. 

Yamamoto, Tanaka, and Murayama 21 provided a detailed examination of the separation bubble at a tail 
incidence of 0 degrees. They showed a large increase in separation from 3 to 4 degrees angle of attack for 
the UPACS structured and TAS unstructured codes. Increasing the chordwise and spanwise grid resolution 
parametrically in the wing root juncture region increased the extent of separation. Using the QCR 22 nonlinear 
eddy viscosity Reynolds stress term eliminated the bubble that was present with a linear eddy viscosity 
Reynolds stress term. In a second paper, Yamamoto, Tanaka, and Murayama 23 examined the CRM with 
and without QCR and found that QCR eliminated the bubble at 4 degrees angle of attack for both the SST 
and Spalart-Allmaras 24 (S-A) turbulence models. Hashimoto et al. 25 applied three different grid topologies 
at 4 degrees angle of attack and concluded, “The separation lines and Cp distribution are largely affected by 
grid topologies.” Refining the region between the wing and tail did not induce a significant change on wing 
root juncture separation, as reported by Hashimoto et al. 26 

Brodersen et al. 2 ' predicted more separation on a prismatic boundary layer grid having a refined wing 
root juncture region than a hexahedral boundary layer grid with a coarser wing root juncture. Crippa 28 
shows both attached and separated wing root juncture flow. The Solar grid has a very coarse wing root 
juncture grid due to a reduction of the boundary layer hexahedra height in corners as a result of the Solar 
generation technique. The SolarChimera5 grid has an extremely well refined wing root juncture grid with 
an orthogonal overset C topology. The SolarChimera5 grid showed wing root separation, but the original 
Solar grid did not show any separation in this region. The SolarChimera5 separation bubble was small 
at the design condition, but grew dramatically in size at 4 degrees angle of attack for the S-A turbulence 
model. The S-A, S-A with rotation/curvature correction, and Reynolds stress turbulence models all had 
the same behavior at the design condition. Crippa also computed a local error estimate of dissipation on 
drag coefficient at the design condition, which indicated that dissipation in the coarse wing root grid of the 
Solar grid had a significant impact on drag. 

Lee-Rausch et al. 29 utilized a fine mixed-element wing root juncture grid and FUN3D with the Venkatakr- 
ishnan limiter. A small wing root separation bubble was shown at the design condition that grew larger with 
grid refinement. Mavriplis and Long 30 utilized the same mixed-element grids with NSU3D. They showed a 
small separation bubble in the wing root juncture region and flow separation behind a strong aft leg of a 
upper surface lambda shock at 4 degrees angle of attack. NSU3D utilizes an edge-based viscous operator, 
which results in a TLNS approximation on orthogonal grids. The Edge solver 31,32 was applied to various 
grid resolutions in wing root juncture and different turbulence models, but no separation was observed. Edge 
also utilizes an edge-based viscous operator. 

Vos, Sachi, and Gehri 33 used 6 different structured multiblock grid systems to compare near-field and 
far-held drag methods. They indicate that separation was predicted with some of the grids, but details of 
the separation region extent were not provided. They indicate, “Grid resolution is an important parameter 
which influence the resolution of these how separations. . . ” Other researchers 34-36 also examined the DPW- 
IV CRM cases, but did not provide details of the wing root juncture how. 

This large body of work performed in conjunction with DPW-IV indicates that the CRM side-of-body 
separation is very sensitive to the CFD scheme used, grid topology, and grid resolution. This sensitivity 
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is especially high in the wing root region. This side-of-body flow separation on the CRM appears to be 
the combination of both modeling and discretization issues. It is difficult to isolate all of the individual 
factors that result in the prediction of separated flow, but some trends appear evident from this review of 
past studies. The modeling options of TLNS and QCR appear to reduce the size of the separation bubble 
as compared to FNS with linear eddy viscosity Reynolds stress term. The use of coarser grids or more 
dissipative convective schemes reduce the separation bubble extent or eliminate this separation completely. 

IV. DPW-V 

The CRM wing-body configuration was examined in the subsequent workshop, DPW-V. DPW-IV identi- 
fied that “Generating a consistent set of grids for the purpose of grid-convergence studies remains a challenge, 
especially for unstructured meshes.” 1 In DPW-V, a common grid system with consistent grid refinement 
factors was required for all participants. 37 The base of the common grid system is a point matched structured 
grid, which was converted for use with overset and unstructured grid flow solvers. This allows a large ma- 
jority of methods used in previous workshops to be applied to consistently refined grids with identical point 
distributions. The variation of computed lift and pitching moment between flow solvers grew dramatically 
between 3 and 4 degrees angle of attack for DPW-IV. 1 To better characterize this flow regime, the angle of 
attack resolution is increased to a quarter of a degree between 2.5 and 4 degrees angle of attack for DPW-V. 
In support of the Turbulence Model Benchmarking Working Group, optional test cases from the Turbulence 
Modeling Resource website 38 were requested. The cases requested for DPW-V are listed here: 

Case 1 — Common Grid Study The 0.5 lift coefficient CRM design condition is simulated on a series of 
DPW committee provided uniformly refined grids to examine grid convergence. 

Case 2 — Buffet Study Cases are requested on the medium DPW committee provided grid at every quar- 
ter of a degree between 2.5 and 4 degrees angle of attack where separation was observed in DPW-IV. 

Case 3 (Optional) — Turbulence Model Verification Three optional cases are specified to study the 
implementation of turbulence models in a controlled study. 

A series of point-matched structured-block grids are provided by the DPW-V committee. 3 ' The com- 
mittee converted the structured grids into an unstructured hexaheclral grid format and formed hybrid grids 
by dividing these hexahedral elements into prisms and tetrahedra. The number of elements in the LI to L5 
grids is shown in Table 2. The L6 grid was not provided by the committee in hybrid form. 


Table 2. Committee supplied structured and hybrid grids. 


Grid Level Nodes 

Hexes 

Prisms 

Tetrahedra 

LI 

660,177 

638,976 

425,984 

2,555,904 

L2 

2,204,089 

2,156,544 

1,437,696 

8,626,176 

L3 

5,196,193 

5,111,808 

3,301,376 

20,766,720 

L4 

17,441,905 

17,252,352 

11,261,952 

69,728,256 

L5 

41,231,169 

40,894,464 

26,411,008 

166,133,760 


V. Method Description 

Three RANS codes are used in this study. CFL3D, FUN3D, and NSU3D have applicability to a large 
range of configurations and provide a large range of options. Only the options used on the CRM in this 
study are described. For the constant coefficient of lift CL grid convergence study, all three codes employ a 
CL driver that allows convergence to a fixed lift coefficient value. Based on user defined inputs, the codes 
periodically adjust the angle of attack during the computation based on the current CL value until the 
prescribed CL value is achieved at convergence. 
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V.A. CFL3D Flow Solver 


CFL3D 39 is a structured-grid multi-zone cell-centered finite-volume method. It was applied to previous 
workshops, including DPW-IV. 19 It uses third-order upwind-biased spatial differencing on the convective 
and pressure terms, and second-order differencing on the viscous terms; it is globally second-order accurate. 
Roe’s flux difference-splitting method 40 is used to obtain fluxes at the cell faces. The option in CFL3D 
to model the FNS mean-flow equations is exercised for all cases. Only the S-A one-equation turbulence 
model was used for this study. The turbulence model diffusion terms use the thin-layer approximation. 
The mean-flow equations are advanced in time with the eddy-viscosity fixed then the turbulence model 
is advanced in time with the mean-flow solution fixed. This alternating turbulence model and mean-flow 
iteration technique is named a loosely coupled iteration. The solution is advanced in time with an implicit 
approximate factorization method. CFL3D employs local time-step scaling, grid sequencing, and multigrid 
to accelerate convergence to steady state. 

V.B. FUN3D Flow Solver 

FUN3D 41,42 is a finite-volume RANS solver in which the flow variables are stored at the vertices or nodes 
of the mesh. FUN3D was used in previous workshops, including DPW-IV. 29 FUN3D solves the equations 
on mixed element grids, including tetrahedra, pyramids, prisms, and hexahedra. At interfaces delimiting 
neighboring control volumes, the inviscicl fluxes are computed using an approximate Riemann solver based 
on the values on either side of the interface. Roe’s flux difference splitting 40 is used in the current study. For 
second-order accuracy, interface values are obtained by a U-MUSCL scheme 43,44 with gradients computed at 
the mesh vertices using an unweighted least-squares technique. The U-MUSCL scheme coefficient x is set to 
0.0 for purely tetrahedral grids and 0.5 for grids with mixed element types. The Venkatakrishnan limiter 45 
is used for some cases in this study. This dimensional limiter is scaled to the mean aerodynamic chord to 
have the same behavior as the airfoil example with unit chord in Venkatakrishnan. 45 

For tetrahedral meshes, the full viscous fluxes are discretized using a finite-volume formulation in which 
the required velocity gradients on the dual faces are computed using the Green-Gauss theorem. On tetra- 
hedral meshes this is equivalent to a Galerkin type approximation. For non-tetrahedral meshes, the same 
Green-Gauss approach can lead to odd-even decoupling. A pure edge-based approach can be used to cir- 
cumvent the odd-even decoupling issue but yields only approximate viscous terms. For non-tetrahedral 
meshes, the edge-based gradients are combined with Green-Gauss gradients; this improves the h-ellipticity 
of the operator and allows the complete viscous stresses to be evaluated. 41,46 This formulation results in 
a discretization of the FNS equations. The option of only using edge-based viscous terms (equivalent to a 
TLNS approximation for orthogonal grids) is also available and is used for a data set in this study. The 
diffusion term in the turbulence model is handled in the same fashion as the mean-flow viscous operator. 

The S-A model may be solved loosely coupled to the mean-flow equations in the same fashion as CFL3D. 
FUN3D also provides a tightly coupled iteration where the turbulence model and mean-flow equations are 
updated simultaneously with a linearization of all coupling terms. The cases for DPW-V were computed 
assuming fully turbulent flow. The negative S-A model 4 ' is available as well as the original model. The 
negative S-A model exhibits better iterative convergence properties, especially in the highly clustered wake 
region of the DPW-V committee provided grids. 

The solution at each time-step is updated with a backwards Euler time-differencing scheme. At each time 
step, the linear system of equations is approximately solved with either a multi-color point-implicit procedure 
or an implicit-line relaxation scheme. This implicit relaxation scheme can be used as a preconditioner 
to Generalized Conjugate Residual (GCR) 48 for stabilization and convergence acceleration. Local time- 
step scaling is employed to accelerate convergence to steady-state. Two other options are exercised on 
some cases to improve iterative convergence rate and robustness: a Jacobian- free Newton method and 
adaptive CFL specification. The Newton method utilized full linearization of the nonlinear residual with a 
pseudo-time addition. It is computed with a first-order Frechet derivative to enhance the standard defect 
correction approach. 49 The solution of this exact linearization is also wrapped with GCR for stabilization and 
convergence acceleration. The CFL specification for the nonlinear psudo-time term and the preconditioner 
can be adjusted with an adaptive approach. 49 
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V.C. NSU3D Flow Solver 


The NSU3D code is an unstructured mesh multigrid RANS solver for high-Reynolds number external aerody- 
namic applications. It has been used in previous workshops, including DPW-IV. 30 The NSU3D discretization 
employs a vertex-based approach where the unknown fluid and turbulence variables are stored at the ver- 
tices of the mesh. The fluxes are computed on faces delimiting dual control volumes with each dual face 
being associated with a mesh edge. This discretization operates on hybrid mixed-element meshes. Generally, 
prismatic elements are employed in highly stretched boundary layer regions and tetrahedral elements are 
employed in isotropic regions of the mesh away from the aircraft surfaces. 

The convective terms are discretized as central differences with added matrix dissipation. Second-order 
accuracy is achieved by formulating these dissipative terms as an undivided biharmonic operator, which is 
constructed in two passes of a nearest-neighbor Laplacian operator. In the matrix form, this dissipation 
is similar to that produced by a Riemann solver gradient-based reconstruction scheme. It is obtained by 
replacing the difference in the reconstructed states on each side of the control volume interface by the undi- 
vided differences along mesh edges resulting from the biharmonic operator construction. These differences 
are then multiplied by the characteristic matrix to obtain the final dissipation terms. 

The baseline NSU3D discretization employs a finite-difference scheme to approximate the thin-layer form 
of the viscous terms for the Navier-Stokes (N-S) equations, although this is done in a multidimensional 
fashion, by computing a Laplacian of the velocity field. 50 The main approximation in this approach is the 
omission of the cross-derivative viscous terms, and the assumption of a locally constant viscosity. NSU3D 
incorporates the S-A turbulence model and the cases for DPW-V were run fully turbulent. 

The basic time-stepping scheme in NSU3D consists of a three stage explicit scheme. Convergence is 
accelerated by a local block- Jacobi preconditioner in regions of isotropic grid cells. In boundary layer regions, 
where the grid is highly stretched, a line preconditioner is employed to relieve the stiffness associated with 
the mesh anisotropy. 51 An agglomeration multigrid algorithm is used to further enhance convergence to 
steady-state. 50,52 The multigrid solver is implemented as a non-linear full approximation scheme solver. 
The Jacobi and line preconditioners are used to drive the various levels of the multigrid sequence, resulting 
in a rapidly converging solution technique. 

At DPW-V, it was noted that NSU3D predicted lower skin friction values than the other methods 
presented. This observation lead to an examination of the skin friction calculation method in NSU3D. 
The grid lines in the DPW-V committee provided grids were determined to be less orthogonal to the solid 
surfaces than typically experienced. To improve the NSU3D DPW-V results, a more general skin-friction 
calculation was implemented in NSU3D which accommodates non-orthogonal grid line intersections with 
viscous surfaces. The updated skin-friction calculation results are provided in this study. 

VI. Test Cases and Results 

The required test cases for DPW-V are a uniform grid convergence study at constant lift with committee 
provided grids and a buffet study on the medium L3 grid between 2.5 and 4 degrees angle of attack. The 
Mach number is 0.85 and the Reynolds number is 5 million per mean aerodynamic chord or 18,129 per full- 
scale inch for all cases. The uniform refinement study is performed for CFL3D, NSU3D, and FUN3D with 
and without the Venkatakrishnan limiter. The uniform refinement results provide an excellent verification 
data set for the FUN3D output-adaptive scheme applied to the LI hybrid grid at the design conditions. 
This allows the adaptive scheme to be compared to a consistently refined set of grids. The buffet study is 
partitioned into applying different CFD schemes to DPW-V L3 medium grids and applying FUN3D to grids 
with different grid topologies. 

VI. A. Grid Convergence Study 

A grid convergence study was performed at the design condition of 0.5 coefficient of lift. The L1-L6 structured 
grids were provided by the DPW committee, but CFL3D failed to reach satisfactory iterative convergence 
on the L6 grid. So, only the L1-L5 results are presented for CFL3D. NSU3D results are only available for 
the L1-L4 grids and use the improved skin friction calculation for non-orthogonal grids. Only L1-L5 grids 
are analyzed with FUN 3D, because the L6 grid was not provided by the committee in hybrid form. 

Drag coefficient CD_TOT is plotted as a function of a characteristic length squared H2 in Fig. 3(a) and 
pitching moment coefficient CM_TOT is plotted as a function of H2 in Fig. 3(b). H2 is computed as the 
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(number of control volumes) -2 / 3 . This exponent is based on two assumptions; the characteristic length of the 
grid varies with the cube root of the cell volume and the error in the solution decreases asymptotically with 
the characteristic length squared for a 2 nd -order method. When these assumptions are met, the computed 
outputs should vary linearly with H2. The legend description FUN3D is FNS without a limiter and FUN3D- 
V is FNS with the Venkatakrishnan limiter. On the finer grids, the difference between FUN3D-V and FUN3D 
diminishes. The CD_TOT values of the finest grid FUN3D, FUN3D-V, and CFL3D are within one count, 
0.0001, of drag coefficient. The finest grid NSU3D simulation (L4) predicted a slightly lower CD_TOT than 
the L4 grid simulations of the other methods. All four methods indicated CM_TOT became more negative 
as the grid was refined. The angle of attack ALPHA required for a coefficient of lift of 0.5 is shown in Fig. 4. 
This ALPHA decreased consistently across the methods as the grid was refined. This decrease in ALPHA 
with grid refinement at constant coefficient of lift is analogous to an increase in coefficient of lift with grid 
refinement at constant ALPHA. 




(a) Coefficient of drag. (b) Coefficient of pitching moment. 

Figure 3. Coefficient of drag and pitching moment as a function of characteristic grid spacing at CL — 0.5 for uniformly 
refined grids. 



Figure 4. Angle of attack as a function of characteristic grid spacing at CL — 0.5 for uniformly refined grids. 


The pressure, CD_PR, and skin friction, CD_SF, components of CD_TOT are shown in Fig. 5. CD_PR 
exhibits a larger variation than CD_SF. The CDJPR is very similar for FUN3D and FUN3D-V on the two 
finest grids. CFL3D has a lower CD_PR than FUN3D and FUN3D-V for all grids, but the difference between 
the methods decreases with grid refinement. The CFL3D CD_SF is extremely insensitive to grid resolution. 
The three unstructured grid methods show a very similar trend of increasing CD_SF with grid refinement. 
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(a) Pressure contributions to coefficient of drag. (b) Skin friction contributions to coefficient of drag. 

Figure 5. Components of coefficient of drag as a function of characteristic grid spacing at CL — 0.5 for uniformly refined 
grids. 


VI. B. Output-Based Grid Adaptation 

Output-based adaptation has been applied to the Shock-Boundary Layer Interaction Workshop 53 and the 
High Lift Prediction Workshop. 54 These previous applications had significant off-body features that allowed 
the use of adaptive mechanics with a frozen boundary layer region. Drag calculation on the CRM configu- 
ration is very sensitive to surface and boundary layer resolution, so a frozen boundary layer approach is not 
suitable. A boundary element subdivision technique was implemented for this study to refine surface and 
boundary layer grids in conjunction with the off-body tetrahedral regions. 

An example of an adapted prismatic and tetrahedral flat plate boundary layer grid is shown in Fig. 6. 
Figure 6(a) shows the original grid and Fig. 6(b) shows an example of an adapted grid to illustrate the 
method. The lower near-body portion of the grid is constructed of stacks of prismatic elements and the 
upper portion outside of the boundary layer prismatic grid contains tetrahedral elements. In Fig. 6(b), the 
prismatic stacks have been subdivided tangentially to the lower surface. 




(a) Original grid. 


(b) Tetrahedral adaptation with prism subdivision. 


Figure 6. Flat plate grid illustrating adaptation technique. 


The prisms stacks are not refined in the normal direction in this application because the initial CRM grid 
met the y + < 1 turbulence model recommendation for a majority of the solid boundaries. A fixed fraction of 
prism element sides with the largest indicator are flagged for refinement. For this flat plate illustration, the 
jump in density over the edge is the indicator. For the CRM, the scalar Venditti 55 remaining error estimate 
is the indicator. The refinement flags are propagated up and down the sides of the prism stacks so that 
all prisms in a stack are refined the same manner. New nodes introduced on the boundary are placed on a 
linear, faceted representation of the surface; the new nodes are not projected to the actual geometry. 

An anisotropic metric specifies the request for the adapted tetrahedra density and orientation. In the 
flat plate illustration, the Mach Hessian scaled by density jumps of Bibb et al. 56 is the specified metric. 
For the CRM, the Venditti 55 output-based metric is specified. The tetrahedral elements are adapted with a 
general metric-based edge primitive scheme. 57,58 Work is ongoing to improve the parallelization of this edge 
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primitive scheme. However, in its current form, the partition boundaries are not evenly adapted. For the 
parallel CRM example, an existing parallel metric-based scheme with element splitting, element collapse, 
element swap, and node movement 59, 60 is also utilized to modify the tetrahedral elements because it has a 
more mature parallelization scheme. 

Figure 7 shows the effects of uniform refinement, FUN3D-V from Fig. 3(a), and drag adaptation on 
CD_TOT. FUN3D-V-METRIC is the existing metric based scheme with frozen prismatic elements. FUN3D- 
V-SUBDIV is the new edge primitive scheme that includes fixed fraction prismatic element subdivision. The 
first four FUN3D-V-METRIC adaptations (without FUN3D-V-SUBDIV) appears to predict a higher drag 
value than uniform refinement, which includes refinement in the boundary layer. The next adaptation 
includes subdivision of the boundary layer prisms (FUN3D-V-SUBDIV) and produces a large reduction 
in drag. The three final FUN3D-V-METRIC adaptations appear to approach the uniformly refined grid 
CD_TOT but do so with a smaller grid. The choice of the number of FUN3D-V-METRIC and FUN3D-V- 
SUBDIV adaptations are arbitrary in this application. Work is ongoing to better integrate the subdivision 
with the metric-based adaptation to retain the benefits of both methods in a more efficient and automated 
manner. 

Lee and Pulliam 61 applied an overset adaptation scheme to the DPW-IV CRM configuration with hor- 
izontal tail, but their adaptation scheme appears to asymptote to larger coefficient of drag values than 
uniformly refined grids. Element subdivision (triggered by the undivided gradient of density) has been used 
on a transport configuration by Mavriplis, 62 but the effect on computed forces was not provided. 

The original and adapted surface grids are shown in Fig. 8. A majority of the wing is refined with the 
nose, tail, and fuselage over the wing. The original and adapted symmetry plane are shown in Figure 9; the 
artifacts of the structured grid generation are shown as a tightly packed grid emanating from the top and 
bottom of the nose and tail. These structured-grid artifacts are removed by the grid adaptation and the 
wake behind the fuselage is refined. The cause of the adapted grid high density clumps is unknown but may 
be due to noise in the current grid spacing estimate or a deficiency of the grid adaptation mechanics. Mach 
number and continuity equation adjoint are shown on the original LI grid in Fig. 10. The grid is refined in 
regions where the adjoint solution indicates drag sensitivity. 



Figure 7. Coefficient of drag as a function of characteristic grid spacing at CL = 0.5 for uniformly refined and adapted 
grids. 


VI. C. Buffet Study (Committee Provided Grids) 

The lift CL_TOT and drag CD_TOT coefficients for the L3 DPW-V committee provided grids are plotted 
as angle of attack ALPHA increases in 0.25 degree increments in Fig. 11. FUN3D-V denotes FNS with 
the Venkatakrishnan limiter, FUN3D denotes FNS with unlimited reconstruction, and FUN3D-TL denotes 
edge-based N-S terms with unlimited reconstruction. CFL3D, NSU3D, and FUN3D-TL show very similar 
trends. FUN3D and FUN3D-V both show dramatic reductions in CL_TOT and CD_TOT, but the limiter 
delays and reduces the magnitude of the reduction. CL MOT is plotted as a function of pitching moment 
coefficient CM_TOT in Fig. 12 where the CM_TOT axis is reversed. CM_TOT trends changes abruptly at 
the same ALPHA as the abrupt CL_TOT change for FUN3D and FUN3D-V. Osusky et al. 63 examined the 
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(a) Original grid. 


(b) Surface grid subdivision. 


Figure 8. CRM surface grid. 



(a) Original grid. 



(b) Metric-based adaptation. 


Figure 9. CRM symmetry plane grid. 



(a) Mach number. 


(b) Continuity equation adjoint solution. 


Figure 10. Mach number and adjoint solution on the CRM symmetry plane. 
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DPW-V L3 committee grid and showed a dramatic drop in lift, drag, and negative pitching moment at 3.25 
degrees angle of attack and above. 



(a) Coefficient of lift as a function of angle of attack. 

Figure 11. 



(b) Coefficient of drag as a function of angle of attack. 


Buffet study total forces for L3 grid. 



Figure 12. Coefficient of lift as a function of pitching moment coefficient for L3 grid. 


The abrupt reductions in forces are accompanied by a significant increase in wing root separation bubble 
size, Fig. 13. Figure 13(a) is the butt line of the spanwise extent of the separation bubble, BLJ3UB, as a 
function of ALPHA. The wing root junction is at butt line 120.267 indicating the bubble is small at low 
ALPHA. Figure 13(b) is the fuselage station of the forward extent of the separation bubble, FSJ3UB. The 
trailing edge of the wing root junction is at the 1458.68 fuselage station and the vertical axis has been 
reversed. A lower fuselage station implies a larger forward bubble extent. 

FUN3D-TL has a very small bubble that remained the same size over the range of ALPHA examined. 
No bubble is reported for CFL3D or NSU3D. The FUN3D and FUN3D-V predicted separation bubbles grow 
with increasing ALPHA. FUN3D shows a dramatic increase in lateral and forward extent at 3.25 degrees 
ALPHA and FUN3D-V shows a similar increase at 3.75 degrees ALPHA. The ALPHA trends of the forces, 
Fig. 11, and moment, Fig. 12, also change at the same ALPHA with rapid separation extent increases. 

The value of the <f> function of the Venkatakrishnan limiter is examined for FUN3D-V to illustrate potential 
differences between the reconstruction of FUN3D and FUN3D-V. Fig. 14 shows the limiter function with 
constant butt line 151.074 and 697.333 volume slices near 13% and 60% span. Fig. 14(a) has two levels of 
nested zoom to show details of the outer edge of the boundary layer and control volume adjacent to the 
solid wall. The upper wing shock, wake, outer boundary layer, and first control volume on the wall are 
areas where the Venkatakrishnan limiter is active (with a value of <j> < 1.0). Although the limiter function 
c t> is close to one in the majority of the field at an angle of attack of 3.25, it has a dramatic effect on the 
side-of-body separation, Fig. 13, forces, Fig. 11, and pitching moment, Fig. 12. 
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(a) Lateral extent of the wing root separation. (b) Longitudinal extent of the wing root separation. 

Figure 13. Wing root separation bubble extent for the buffet study on the DPW committee grids. No separation 
reported for CFL3D or NSU3D. 




Figure 14. 


(a) 77 = 0.1306. (b) 77 = 0.6028. 

The value of the Venkatakrishnan limiter, </>, for the x-momentum equation at a. = 3.25. 
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Coefficient of skin friction color contours and coefficient of pressure contour lines in 0.1 increments are 
shown in Fig. 15. The series of subplots in Fig. 15 detail the separation bubble growth for FUN3D FNS with 
the unlimited scheme on the L3 hybrid grid. At ALPHA= 3.25 and lower, a lambda shock pattern is seen on 
the upper surface of the wing in the coefficient of pressure contour lines. The skin friction is low in the wing 
root junction behind the aft leg of the lambda shock. The wing root separation is contained in this low skin 
friction region. The skin friction has also decreased behind the outer wing shock. At ALPHA= 3.50 and 
above, the wing root separation region has moved to the forward leg of the lambda shock. The displacement 
of the streamlines due to the separation region has strengthened the forward leg of the lambda shock (the 
coefficient of pressure contour lines are closer together) and nearly eliminated the aft leg. The outer wing skin 
friction increases with the change in upper wing shock topology. The skin friction then continues to decreases 
behind the outer wing shock as angle of attack is increased. The upper wing shock locations (clusters of 
coefficient of pressure contour lines) appear to correlate with the initiation of wing root separation. 


VI. D. Buffet Study (Custom Grids) 

In addition to the L3 DPW-V grid, grids based on a DPW-IV grid and custom grids are analyzed for the 
buffet study. The FNS unlimited reconstruction FUN3D scheme is used on all of the grids in this section. 
This allows these different grid topologies to be evaluated with the same scheme. The number of nodes 
and elements in these custom grids are listed in Table 3. Both tetrahedral grids, DPW4T and CUSTT, 
are generated with VGRID. 64 The mixed-element grids, DPW4M and CUSTM, are created by merging 
the boundary layer tetrahedral elements. 52 A small number of Steiner nodes are required to complete the 
tetrahedra merging process, 52 so the mixed-element grids have slightly more nodes than the orignal purely 
tetrahedral grids. 

The DPW4T tetrahedral and DPW4M mixed-element grids are based on the the mixed-element medium 
no tail grid described by Lee-Rausch et al. 29 where the FUN3D calculations utilized FNS and the Venkatakr- 
ishnan limiter. The VGRID sources associated with the horizontal tail have been omitted to coarsen the 
grid in the tail region but retain the wing grid resolution of Lee-Rausch et al. 29 The CUSTT and CUSTM 
grids have less resolution in the wing root junction, but are larger overall with less spanwise stretching. The 
FUN3D-CUSTT and FUN3D-CUSTM custom grids have twice the total node count of the DPW4T and 
DPW4M grids and three times the total node count of the DPW-V L3 grid. 


Table 3. Other grids examined. 


Nodes Prisms 


Grid 

DPW4T 

DPW4M 

CUSTT 

CUSTM 


7,036,245 

7,053,021 

14,579,511 

14,625,035 


0 

10,520,438 

0 

13,321,168 


Pyramids Tetrahedra 
0 41,497,260 

59,938 9,899,950 

0 85,524,780 

102,366 45,584,164 


The buffet study was repeated with these custom grids. CL_TOT and CD_TOT are plotted as ALPHA 
is increased in 0.25 degree increments in Fig. 16. CL_TOT is plotted as a function of CM_TOT in Fig. 17. 
The CM_TOT axis is reversed. The unlimited reconstruction FNS FUN3D calculation on the L3 hybrid grid 
is included for comparison. The mixed-element FUN3D-CUSTM produced larger CL_TOT and CD_TOT 
and a more negative CM_TOT than the tetrahedral FUN3D-CUSTT for the entire ALPHA range examined. 
FUN3D-DPW4T shows a large decrease in lift at 3 degrees angle of attack. FUN3D-DPW4M shows a large 
decrease in lift at 2.75 degrees angle of attack. 

As with the workshop provided grids, the abrupt reductions in forces are accompanied by a significant 
increase in wing root separation bubble size, Fig. 18. Figure 18(a) is the butt line of the spanwise extent of 
the separation bubble, BLJ3UB, as a function of ALPHA. The wing root junction is at a butt line of 120.267 
indicating the bubble is small at low ALPHA. Figure 18(b) is the fuselage station of the forward extent of 
the separation bubble, FSJ3UB. The trailing edge of the wing root junction is at the 1458.68 fuselage station 
and the vertical axis has been reversed. A lower fuselage station implies a larger forward bubble extent. 

FUN3D-CUSTT has a very small bubble that remained the same size over the range of ALPHA examined. 
FUN3D-CUSTM also has a small bubble with a small increase in size with angle of attack. FUN3D, FUN3D- 
DPW4T, and FUN3D-DPW3M predicted separation bubbles that grow with increasing ALPHA. FUN3D- 
DPW4T and FUN3D-DPW4M have the largest separation bubble extents for the entire ALPHA range, which 
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(a) ALPHA= 2.75. 


(b) ALPHA= 3.00. 




(c) ALPHA= 3.25. 


(d) ALPHA= 3.50. 




(e) ALPHA= 3.75. 


(f) ALPHA= 4.00. 


Figure 15. X-component of skin friction color contours and coefficient of pressure contour lines for FUN3D FNS with 
the unlimited scheme on the L3 hybrid grid. 
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ALPHA 


ALPHA 


(a) Coefficient of lift as a function of angle of attack. 


(b) Coefficient of drag as a function of angle of attack. 


Figure 16. Buffet study total forces for L3 and custom grids. 



Figure 17. Coefficient of lift as a function of pitching moment coefficient for L3 and custom grids. 
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starts near the design condition of CL_TOT= 0.5. A dramatic increase in lateral and forward extent is shown 
by FUN3D at 3.25 degrees ALPHA, FUN3D-DPW4T at 3.00 degrees ALPHA, and FUN3D-DPW4M at 2.75 
degrees ALPHA. The ALPHA trends of the forces, Fig. 16, and moment, Fig. 17, also change at the same 
ALPHA with rapid separation extent increases. 



(a) Lateral extent of the wing root separation. 



(b) Longitudinal extent of the wing root separation. 


Figure 18. Wing root separation bubble extent for the buffet study for L3 and custom grids. 


The inboard portion of the upper wing is shown in Fig. 19. The surface grids and streamlines for the 
mixed element grids at ALPHA= 4.00 are shown for the unlimited reconstruction FNS FUN3D. The grid 
resolution is fine at the leading and trailing edges of all three grids. The DPW4T and DPW4M grids, 
Fig. 19(c), have the finest wing root grid (the lower portion of the subfigure). The L3, Fig. 19(a), grid has 
slightly less wing root resolution. The CUSTT and CUSTM grids, Fig. 19(e), has the coarsest wing root 
grid. The grids with more wing root grid resolution were also more likely to have wing root separation at 
lower angles of attack in Fig. 18. 

The large FUN3D and FUN3D-DPW3M separation bubbles traced by the streamline have similar struc- 
ture and extent in Fig. 19(b) and Fig. 19(d). The stronger shock locations can be observed initiating the 
abrupt outward turning of the streamlines. FUN3D-CUSTM has a small separation bubble isolated to the 
training edge near the fuselage in Fig. 19(f). The strong shock in the aft leg of the FUN3D-CUSTM upper 
surface lambda shock induces an abrupt outward turning of the streamlines, which is very clear near the 
trailing edge Yehudi break. Osusky et al. 63 examined the L3 grid and showed streamlines with a large 
separation bubble of similar topology at 3.25 degrees angle of attack and above. 

VII. Buffet Study Wing Root Separation Bubble 

All the methods showed no separation bubble or a small bubble at the design condition of CL= 0.5. The 
large separation bubbles where observed in the higher angle of attack buffet study. The initiation of the 
separation region appears to be connected to the location of the upper wing shocks. When the separation 
bubble grows dramatically, as angle of attack increases, it tends to move from the aft leg of the upper wing 
lambda shock to the forward leg. Simultaneously, the strongest leg of the lambda shock changes from the 
aft leg to the forward leg. This shock and separation bubble interaction mechanism, as seen in this study 
and a review of DPW-IV publications, appears to be the dominate mechanism. The resultant shock, surface 
pressure, skin friction, and streamline topologies of the no/small separation bubble attractor and the large 
separation bubble attractor appear to be similar across a wide range of discretizations. However, it is less 
obvious how to predict if a particular discretization with exhibit the no, small, or large separation bubble 
topology. 

CFL3D did not show separation with FNS, but has a very different dimension-by-dimension reconstruc- 
tion scheme than the unstructured FUN3D and NSU3D. There is no reported separation bubble for NSU3D, 
which utilizes edge-based viscous terms that produce the TLNS approximation on orthogonal grids. FUN3D 
showed a wide range of separation bubble size for the various discretization options and grid topologies 
examined. 

On the DPW-V L3 grid, FUN3D with unlimited reconstruction and FNS showed a large increase in 
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(c) DPW4T and DPW4M. 


(d) FUN3D-DPW4M. 


(e) CUSTT and CUSTM. 


(f) FUN3D-CUSTM. 


Figure 19. Inboard section of the upper wing surface near the fuselage, surface grid and streamlines at ALPHA- 4.00, 
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separation bubble size and drop in CL at ALPHA= 3.25. FUN3D with the Venkatakrishnan limiter delayed 
the ALPHA with the large separation bubble size growth and reduced its lateral extent. The unlimited 
reconstruction and edge-based viscous operator producing a TLNS approximation resulted in the small 
separation bubble topology at ALPHA= 4.00. 

FUN3D with unlimited reconstruction and FNS was applied to grids with different grid topologies to 
examine the effect of wing root grid refinement and separation bubble topology. A DPW-IV based grid 
with finer wing root grid resolution than the DPW-V L3 grid transitioned to the large separation bubble 
topology at a lower ALPHA than the L3 grid. This DPW-IV based grid also had a larger bubble near the 
design condition of CL= 0.5. A custom grid with less wing root grid resolution then the DPW-V L3 grid 
did not transition to the large separation bubble topology. The mixed-element versions of each of these 
grids had either transitioned at a lower ALPHA or had a larger bubble than the purely tetrahedral versions. 
Yamamoto, Tanaka, and Murayama 21 also showed the trend of separation bubble growth with wing root 
grid refinement and a reduction in dissipation. 

VIII. Conclusions 

The CFL3D, FUN3D, and NSU3D computational methods were applied to the DPW-V CRM configu- 
ration. These included a structured grid method and two node-based mixed-element methods. The uniform 
grid refinement study at the design lift condition showed a tight grouping in predicted drag for the finest 
grids. Drag-based adaptation was applied at the design lift condition. This included a combination of 
tetrahedral metric-based adaptation and mixed-element subdivision. The grid adapted drag results from a 
coarser grid approached those generated using a finer uniformly refined grid. 

The buffet study using the committee provided grids resulted in a larger range of predicted forces and 
pitching moment coefficients than at the design condition. Abrupt reductions in lift and drag coefficients 
was correlated with an increase in side-of-body separation bubble size as angle of attack was increased. The 
initiation of the wing root separation appears to be related to upper wing shock locations. Two distinct 
topologies emerged: a no/small separation bubble and a large separation bubble. 

It does not appear possible to predict, with absolute certainty, the separation bubble topology of a 
particular discretization based on its properties, but the following trends were observed. The use of the 
Venkatakrishnan limiter delayed the onset angle of attack of this effect as well as its severity. The use 
of edge-based viscous terms reduced the size of the bubble and resulted in higher lift and drag values 
without abrupt increases in side-of-body separation extent as compared to modeling FNS. However, the 
TLNS approximation changes modeling assumptions and may not be appropriate as a predictive tool. 

The participant provided grids had a larger range of grid sizes and differences in predicted force and mo- 
ment coefficients. Converting a tetrahedral grid into a mixed element grid increased the predicted lift, drag, 
and separation bubble extent, which may be the result of a reduction in dissipation over the purely tetra- 
hedral grids. Grids with fine resolution of the wing root juncture region showed more wing root separation 
and where more likely to show abrupt increases in separation region extent. 

In this study, the lowest dissipation scheme on the finest wing root grid resolution predicted the large 
separation bubble topology for FNS with the S-A turbulence model. This indicates that the large separa- 
tion bubble topology may be the correct solution to the RANS equations with the S-A turbulence model. 
However, this large separation bubble topology was not observed in wind tunnel tests indicating that there 
is still a need to improve modeling techniques to predict the wind tunnel results for this configuration and 
possibility other wing root juncture flows. 
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